#for flowchart
library(DiagrammeR)
#to export flowchart as .png
#webshot::install_phantomjs()
library(DiagrammeRsvg)
library(rsvg)
library(grid)
library(xtable)
library(dplyr)
library(ggplot2)
library(gridExtra)
#nest/unnest
library(tidyr)
#map function (kind of like a for loop)
library(purrr)
#tidy model summary
library(broom)
library(readxl)
#for tableby
library(arsenal)

Define geom_split_violin()

Based on https://debruine.github.io/posts/plot-comparison/

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, draw_group = function(self, data, ..., draw_quantiles = NULL){
  data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
  grp <- data[1,'group']
  newdata <- plyr::arrange(transform(data, x = if(grp%%2==1) xminv else xmaxv), if(grp%%2==1) y else -y)
  newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
  newdata[c(1,nrow(newdata)-1,nrow(newdata)), 'x'] <- round(newdata[1, 'x']) 
  if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
    stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <= 
                                              1))
    quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
    aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
    aesthetics$alpha <- rep(1, nrow(quantiles))
    both <- cbind(quantiles, aesthetics)
    quantile_grob <- GeomPath$draw_panel(both, ...)
    ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
  }
  else {
    ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
  }
})

geom_split_violin <- function (mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                               draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, show.legend = NA, 
                               inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, position = position, 
        show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Load Data

load('Data/DataWithPropensities_seed1.RData')

#convert PrimaryDiagnosis to factor
dat3$PrimaryDiagnosis <- factor(dat3$PrimaryDiagnosis, levels = c("Autism", "None"))
levels(dat3$PrimaryDiagnosis)
## [1] "Autism" "None"
dat3$PrimaryDiagnosis <- relevel(dat3$PrimaryDiagnosis, "None")
levels(dat3$PrimaryDiagnosis) = c("TD", "ASD")

Reshape data to combine motion QC levels

#create dummy factor to include all subjects
dat3$noExclusion <- ifelse(dat3$ID>0, "Pass", "Pass")
dat3$noExclusion <- factor(dat3$noExclusion)

#convert KKI_criteria to factor level "Pass"
dat3$KKI_criteria <- factor(dat3$KKI_criteria, levels=c("Pass", "Fail"))

#convert Ciric_length to factor with reference level "Pass" to match KKI_criteria
dat3$Ciric_length <- factor(dat3$Ciric_length, levels=c("Pass", "Fail"))

#combine Ciric_length, KKI, and noExclusion exclusion into one variable
allVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "Ciric_length", "KKI_criteria", "noExclusion", 
                 "PANESS.TotalOverflowNotAccountingForAge", "SRS.Score", "WISC.GAI", 
                 "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw",  "ADOS.Comparable.Total", 
                 "CurrentlyOnStimulants", "HeadCoil", "Sex", 
                 "Diagnosis4Levels", "ADHD_Secondary", "SES.Family", 
                 "Race2", "handedness", "CompletePredictorCases", "YearOfScan", "MeanFramewiseDisplacement.KKI")

idVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan",
                "PANESS.TotalOverflowNotAccountingForAge", "SRS.Score", "WISC.GAI", 
                "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw", "ADOS.Comparable.Total", 
                "CurrentlyOnStimulants", "HeadCoil", "Sex", 
                "Diagnosis4Levels", "ADHD_Secondary", "SES.Family", 
                "Race2", "handedness", "CompletePredictorCases", "YearOfScan", "MeanFramewiseDisplacement.KKI")

qcMelt <- reshape2::melt(dat3[, allVariables],
                         id.vars=names(dat3)[which(names(dat3) %in%  idVariables)], 
                         variable.name = "Motion.Exclusion.Level",
                         value.name = "Included")
## Warning: attributes are not identical across measure variables; they will be
## dropped
#rename exclusion levels
#NOTE: need None to be highest level for geom_split_violin
levels(qcMelt$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

#convert Included to factor with pass as reference
qcMelt$Included <- factor(qcMelt$Included, levels=c("Pass", "Fail"))

#rename levels of value
levels(qcMelt$Included) <- c("Included", "Excluded")

with(qcMelt, table(PrimaryDiagnosis, Motion.Exclusion.Level, Included))
## , , Included = Included
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     151     308  372
##              ASD     29     114  173
## 
## , , Included = Excluded
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     221      64    0
##              ASD    144      59    0
qcMelt$Diagnosis4Levels <- factor(qcMelt$Diagnosis4Levels)

Check consistency of diagnosis variables created by Ben and Liwei

#limit to complete predictor cases
dat3 <- filter(dat3, CompletePredictorCases==1)

#convert Diagnosis4Levels to factor
dat3$Diagnosis4Levels <- factor(dat3$Diagnosis4Levels, levels = c("Autism", "Autism+ADHD", "None"))
table(dat3$Diagnosis4Levels)
## 
##      Autism Autism+ADHD        None 
##          49         102         353
table(dat3$ADHD_Secondary)
## 
##   0   1 
## 402 102

49 Autism + 353 None = 402 ADHD_Secondary=0

1. Impact of motion QC on sample size

tabN<- tableby(PrimaryDiagnosis~Motion.Exclusion.Level, data=filter(qcMelt, Included=="Included"))

summary(tabN,  
        title='N',digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)
N
TD (N=831) ASD (N=316)
Motion.Exclusion.Level
   Strict 151 (18.2%) 29 (9.2%)
   Lenient 308 (37.1%) 114 (36.1%)
   None 372 (44.8%) 173 (54.7%)

1b. Summarize sociodemographic info for complete cases for paper

#sex
qcMelt$Sex <- relevel(as.factor(qcMelt$Sex), "M")

tabSex<- tableby( PrimaryDiagnosis~Sex, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

#age
tabAge<- tableby( PrimaryDiagnosis~AgeAtScan, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

tab12 <- merge(tabSex, tabAge)

#Race
tabRace<- tableby( PrimaryDiagnosis~Race2, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="fe", total=FALSE))

tab123 <- merge(tab12, tabRace)

#SES
tabSES<- tableby( PrimaryDiagnosis~SES.Family, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

tab1234 <- merge(tab123, tabSES)

#Handedness
qcMelt$handedness <- factor(qcMelt$handedness)
qcMelt$handedness <- relevel(qcMelt$handedness, "Right")

tabHand<- tableby( PrimaryDiagnosis~handedness, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

tab12345 <- merge(tab1234, tabHand)

#Currently on Stimulants
qcMelt$CurrentlyOnStimulants <- factor(qcMelt$CurrentlyOnStimulants)

tabStim<- tableby( PrimaryDiagnosis~CurrentlyOnStimulants, data=filter(qcMelt, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

tab123456 <- merge(tab12345, tabStim)
summary(tab123456)
TD (N=372) ASD (N=173) p value
Sex < 0.001
   M 258 (69.4%) 148 (85.5%)
   F 114 (30.6%) 25 (14.5%)
AgeAtScan 0.902
   Mean (SD) 10.383 (1.246) 10.382 (1.356)
   Range 8.020 - 12.980 8.010 - 12.990
Race2 0.001
   N-Miss 0 1
   African American 40 (10.8%) 11 (6.4%)
   Asian 28 (7.5%) 3 (1.7%)
   Biracial 47 (12.6%) 14 (8.1%)
   Caucasian 257 (69.1%) 144 (83.7%)
SES.Family < 0.001
   N-Miss 8 9
   Mean (SD) 54.243 (9.371) 51.396 (9.506)
   Range 18.500 - 66.000 27.000 - 66.000
handedness 0.259
   N-Miss 1 0
   Right 333 (89.8%) 147 (85.0%)
   Left 19 (5.1%) 12 (6.9%)
   Mixed 19 (5.1%) 14 (8.1%)
CurrentlyOnStimulants < 0.001
   N-Miss 1 8
   0 371 (100.0%) 105 (63.6%)
   1 0 (0.0%) 60 (36.4%)

Generate version of table to paste into overleaf

tab <- xtable(as.data.frame(summary(tab123456)))
print(tab, type="latex")
## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Mon Sep 27 21:42:05 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllll}
##   \hline
##  &  & TD (N=372) & ASD (N=173) & p value \\ 
##   \hline
## 1 & **Sex** &  &  & $<$ 0.001 \\ 
##   2 & \&nbsp;\&nbsp;\&nbsp;M & 258 (69.4\%) & 148 (85.5\%) &  \\ 
##   3 & \&nbsp;\&nbsp;\&nbsp;F & 114 (30.6\%) & 25 (14.5\%) &  \\ 
##   4 & **AgeAtScan** &  &  & 0.902 \\ 
##   5 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 10.383 (1.246) & 10.382 (1.356) &  \\ 
##   6 & \&nbsp;\&nbsp;\&nbsp;Range & 8.020 - 12.980 & 8.010 - 12.990 &  \\ 
##   7 & **Race2** &  &  & 0.001 \\ 
##   8 & \&nbsp;\&nbsp;\&nbsp;N-Miss & 0 & 1 &  \\ 
##   9 & \&nbsp;\&nbsp;\&nbsp;African American & 40 (10.8\%) & 11 (6.4\%) &  \\ 
##   10 & \&nbsp;\&nbsp;\&nbsp;Asian & 28 (7.5\%) & 3 (1.7\%) &  \\ 
##   11 & \&nbsp;\&nbsp;\&nbsp;Biracial & 47 (12.6\%) & 14 (8.1\%) &  \\ 
##   12 & \&nbsp;\&nbsp;\&nbsp;Caucasian & 257 (69.1\%) & 144 (83.7\%) &  \\ 
##   13 & **SES.Family** &  &  & $<$ 0.001 \\ 
##   14 & \&nbsp;\&nbsp;\&nbsp;N-Miss & 8 & 9 &  \\ 
##   15 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 54.243 (9.371) & 51.396 (9.506) &  \\ 
##   16 & \&nbsp;\&nbsp;\&nbsp;Range & 18.500 - 66.000 & 27.000 - 66.000 &  \\ 
##   17 & **handedness** &  &  & 0.259 \\ 
##   18 & \&nbsp;\&nbsp;\&nbsp;N-Miss & 1 & 0 &  \\ 
##   19 & \&nbsp;\&nbsp;\&nbsp;Right & 333 (89.8\%) & 147 (85.0\%) &  \\ 
##   20 & \&nbsp;\&nbsp;\&nbsp;Left & 19 (5.1\%) & 12 (6.9\%) &  \\ 
##   21 & \&nbsp;\&nbsp;\&nbsp;Mixed & 19 (5.1\%) & 14 (8.1\%) &  \\ 
##   22 & **CurrentlyOnStimulants** &  &  & $<$ 0.001 \\ 
##   23 & \&nbsp;\&nbsp;\&nbsp;N-Miss & 1 & 8 &  \\ 
##   24 & \&nbsp;\&nbsp;\&nbsp;0 & 371 (100.0\%) & 105 (63.6\%) &  \\ 
##   25 & \&nbsp;\&nbsp;\&nbsp;1 & 0 (0.0\%) & 60 (36.4\%) &  \\ 
##    \hline
## \end{tabular}
## \end{table}

1c. Create exclusion flowchart

1d. Proportion excluded

Define theme for proportion excluded plots

#Easier to make three separate panels and combine them than to use faceting function
My_Theme_prop = theme_light()+theme(
  legend.title =element_blank(),
  axis.title.x = element_text(size = 9),
  axis.title.y = element_text(size = 9),
  plot.title = element_text(size = 30),
  axis.text.x = element_text(size = 6),
  axis.text.y = element_text(size = 8),
  strip.text.x = element_text(size = 10,color="black"),
  strip.background = element_rect(fill="white"))

1d.1 Proportion excluded in each Diagnostic Group

levels(qcMelt$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

motion <- filter(qcMelt, Motion.Exclusion.Level!="None")
motion$Motion.Exclusion.Level <- droplevels(motion$Motion.Exclusion.Level)

levels(motion$PrimaryDiagnosis)
## [1] "TD"  "ASD"
motion <- group_by(motion, PrimaryDiagnosis, Motion.Exclusion.Level, Included)

dx_proportions <- ggplot(motion, aes(x = PrimaryDiagnosis, fill = Included))+
  geom_bar(position = "fill", alpha = .6)+
  facet_grid(~Motion.Exclusion.Level)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme_prop+theme(legend.title =element_blank())+
  theme(legend.title = element_blank())+
  ylab("Proportion of Children")+
  theme(legend.position = "right")+
  theme(legend.margin = margin(t=0, r=0, b=-1, l=-1))+
  theme(legend.key.size=unit(.08, "in"))+  
  theme(axis.title.x = element_blank())+
  theme(axis.text.x = element_text(size = 10))
  
png("Figures/fig_propExcludedDxOnly_cc.png",width=3,height=2,units="in",res=200)
dx_proportions
invisible(dev.off())

dx_proportions

dx_proportions <- dx_proportions +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(size = 6), axis.title.x = element_text(size = 9))+
  xlab("Diagnosis Group")

#Pearson's chi squared tests
extib <- tibble(motion)

exNest <- extib %>% 
  select(c("PrimaryDiagnosis",
           "Motion.Exclusion.Level", "Included")) %>% 
  group_by(Motion.Exclusion.Level) %>% 
  tidyr::nest()

#nested models
ex_chisq <- exNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$PrimaryDiagnosis, .x$Included)))) %>% 
  unnest(stats) 

ex_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic  p.value parameter method 
##   <fct>                  <list>                 <dbl>    <dbl>     <int> <chr>  
## 1 Strict                 <tibble [504 × 2]>      19.9  8.07e-6         1 Pearso…
## 2 Lenient                <tibble [504 × 2]>      10.3  1.30e-3         1 Pearso…

A larger proportion of autistic children were excluded compared to typically developing children using both levels of motion QC.

1d.2 Proportion excluded by ADHD comorbidity in the ASD group

comorbid <- filter(motion, PrimaryDiagnosis=="ASD")
comorbid$PrimaryDiagnosis <- droplevels(comorbid$PrimaryDiagnosis)
comorbid$Diagnosis4Levels <- droplevels(comorbid$Diagnosis4Levels)
levels(comorbid$Diagnosis4Levels)
## [1] "Autism"      "Autism+ADHD"
comorbid <- filter(comorbid, !is.na(Diagnosis4Levels))

comorbid_proportions <- ggplot(comorbid, aes(x = Diagnosis4Levels, fill = Included))+
  geom_bar(position = "fill", alpha = .6)+
  facet_grid(~Motion.Exclusion.Level)+
  scale_x_discrete(labels = c("No", "Yes"))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme_prop+theme(legend.title =element_blank())+
  xlab("Comorbid ADHD?")+
  theme(legend.title = element_blank())+
  #theme(legend.position="top")+
  theme(legend.text = element_text(size=7))+
  ylab("Proportion of Autistic Children")+
  theme(legend.margin = margin(t=0, r=0, b=-1, l=-1))+
  theme(legend.key.size=unit(.1, "in"))+
  theme(legend.position = "none")
  #theme(plot.margin = unit(c(1, 6, 1, 6), "cm"))

prop_legend = cowplot::get_legend(comorbid_proportions + guides(color = guide_legend(ncol = 1))+
                                    theme(legend.position = "right", legend.text = element_text(size = 9),
                                          legend.key.size=unit(.1, "in")))

asdNest <- tibble(comorbid) %>% 
  select(c("Diagnosis4Levels", "CurrentlyOnStimulants",
           "Motion.Exclusion.Level", "Included")) %>% 
  reshape2::melt(id.vars = c("Motion.Exclusion.Level", "Included"),
       measure.vars = c("Diagnosis4Levels", "CurrentlyOnStimulants")) %>% 
  group_by(Motion.Exclusion.Level, variable) %>% 
  tidyr::nest()
## Warning: attributes are not identical across measure variables; they will be
## dropped
#nested models
asd_chisq <- asdNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$value, .x$Included)))) %>%
  unnest(stats)

asd_chisq
## # A tibble: 4 × 7
## # Groups:   Motion.Exclusion.Level, variable [4]
##   Motion.Exclusion… variable   data    statistic p.value parameter method       
##   <fct>             <fct>      <list>      <dbl>   <dbl>     <int> <chr>        
## 1 Strict            Diagnosis… <tibbl…     0.231  0.631          1 Pearson's Ch…
## 2 Lenient           Diagnosis… <tibbl…     0.463  0.496          1 Pearson's Ch…
## 3 Strict            Currently… <tibbl…     6.40   0.0114         1 Pearson's Ch…
## 4 Lenient           Currently… <tibbl…     6.39   0.0115         1 Pearson's Ch…

1d.3 Proportion excluded by Currently on stimulants in ASD group

stim <- filter(motion, PrimaryDiagnosis=="ASD")
stim$PrimaryDiagnosis <- droplevels(stim$PrimaryDiagnosis)
stim$CurrentlyOnStimulants <- factor(stim$CurrentlyOnStimulants)
levels(stim$CurrentlyOnStimulants)
## [1] "0" "1"
levels(stim$CurrentlyOnStimulants) = c("No", "Yes")

stim <- filter(stim, !is.na(CurrentlyOnStimulants))

stim_proportions <- ggplot(stim, aes(x = CurrentlyOnStimulants, fill = Included))+
  geom_bar(position = "fill", alpha = .6)+
  facet_grid(~Motion.Exclusion.Level)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme_prop+
  xlab("Currently On Stimulants?")+
  theme(legend.title = element_blank())+
  #theme(legend.position="top")+
  theme(legend.text = element_text(size=7))+
  ylab("Proportion of Autistic Children")+
  theme(legend.margin = margin(t=0, r=0, b=-1, l=-1))+
  theme(legend.key.size=unit(.1, "in"))+
  theme(legend.position = "none")
  #theme(plot.margin = unit(c(1, 6, 1, 6), "cm"))

1d.4 Combine proportion plots into one figure & save

2. rs-fMRI exclusion probability as a function of age and symptom severity

Covariates

phenoVariables <- c("ID", "PrimaryDiagnosis", 
                    "ADOS.Comparable.Total", 
                    "SRS.Score", 
                    "PANESS.TotalOverflowNotAccountingForAge", 
                    "DuPaulHome.InattentionRaw",
                    "DuPaulHome.HyperactivityRaw", 
                    "AgeAtScan", 
                    "WISC.GAI",
                    "Motion.Exclusion.Level", "Included")

phenoIDs <- c("ID", "PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")

aim1 <- reshape2::melt(qcMelt[, phenoVariables],
                         id.vars=names(qcMelt)[which(names(qcMelt) %in% phenoIDs)])

levels(aim1$variable) <- c("ADOS", "SRS", "Motor Overflow", "Inattention", 
                           "Hyperactivity", "Age", "GAI")

aim1G <- group_by(aim1, PrimaryDiagnosis, Motion.Exclusion.Level, Included, variable)

2a. Univarate GAMs

aim1$delta = rep(NA,length=nrow(aim1))
aim1$delta = ifelse(aim1$Included=="Included",1,0)
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level!="None"))
aim1tib$Motion.Exclusion.Level <- droplevels(aim1tib$Motion.Exclusion.Level)

aim1Nest <- aim1tib %>% 
  group_by(Motion.Exclusion.Level, variable) %>% 
  tidyr::nest()

#nested models
nested_gams <- aim1Nest %>% 
  mutate(model = map(data, ~mgcv::gam(1-delta~s(value, k=-1), data = na.omit(.x), 
                                      family=binomial(link=logit), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE),
         Rsq = map_dbl(model, ~summary(.)$r.sq)) %>% 
  unnest(coefs)

#Ben: correct for 7 lenient and 7 strict 
nested_gams_len <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Lenient")

nested_gams_len$p.fdr = p.adjust(nested_gams_len$p.value, method = "BH")

nested_gams_strict <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Strict")
  
nested_gams_strict$p.fdr = p.adjust(nested_gams_strict$p.value, method = "BH")

#combine to print
nested_gams <- rbind(nested_gams_len, nested_gams_strict)

#list adjusted p values
nested_gams[, c(1:2,6:11)]
## # A tibble: 14 × 8
## # Groups:   Motion.Exclusion.Level, variable [14]
##    Motion.Exclusion.… variable     edf ref.df statistic  p.value    Rsq    p.fdr
##    <fct>              <fct>      <dbl>  <dbl>     <dbl>    <dbl>  <dbl>    <dbl>
##  1 Lenient            ADOS        2.04   2.47     18.4   2.08e-4 0.0442  4.86e-4
##  2 Lenient            SRS         1.83   2.29     17.4   2.95e-4 0.0466  5.17e-4
##  3 Lenient            Motor Ove…  1.00   1.00     15.6   7.58e-5 0.0308  4.54e-4
##  4 Lenient            Inattenti…  1.38   1.67      7.44  1.17e-2 0.0157  1.17e-2
##  5 Lenient            Hyperacti…  2.15   2.70     13.2   4.59e-3 0.0242  5.36e-3
##  6 Lenient            Age         1.00   1.00      9.33  2.26e-3 0.0164  3.17e-3
##  7 Lenient            GAI         1.00   1.00     14.7   1.30e-4 0.0288  4.54e-4
##  8 Strict             ADOS        1.00   1.00     21.9   2.79e-6 0.0444  9.76e-6
##  9 Strict             SRS         1.00   1.00     22.7   1.53e-6 0.0567  9.76e-6
## 10 Strict             Motor Ove…  1.00   1.00     10.2   1.42e-3 0.0194  1.99e-3
## 11 Strict             Inattenti…  1.00   1.00     17.8   2.47e-5 0.0360  4.31e-5
## 12 Strict             Hyperacti…  1.88   2.35     25.0   1.30e-5 0.0516  3.04e-5
## 13 Strict             Age         1.76   2.20      7.73  2.84e-2 0.0129  2.84e-2
## 14 Strict             GAI         1.00   1.00      7.55  6.02e-3 0.0130  7.02e-3
nested_gams <- nested_gams %>% 
           mutate(LB = map(data, ~round(min(na.omit(.x$value)))),
                  UB = map(data, ~round(max(na.omit(.x$value)))),
                  range = map2(LB, UB, ~seq(from=.x, to=.y, by=1)),
                  logpredict = map2(model, range, ~predict(.x, newdata = data.frame(value = .y), type="link",se.fit=TRUE)),
                  fit = map(logpredict, ~plogis(.x$fit)),
                  lCI = map(logpredict, ~plogis(.x$fit-1.96*.x$se.fit)),
                  hCI = map(logpredict, ~plogis(.x$fit+1.96*.x$se.fit)))

Define theme for gam plots

2a.1 Probability of exclusion as a function of ADOS

ados <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_ados <- ggplot(ados, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='Probability of Exclusion', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))

2a.2 Probability of exclusion as a function of SRS

srs <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_srs <- ggplot(srs, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+
  gam_theme+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.3 Probability of exclusion as a function of Inattention

ina <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_in <- ggplot(ina, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.4 Probability of exclusion as a function of Hyperactivity

hi <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_hi <- ggplot(hi, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Hyperactivity")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.5 Probability of exclusion as a function of Motor Overflow

mo <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_mo <- ggplot(mo, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.6 Probability of exclusion as a function of Age

age<- nested_gams %>% 
  filter(variable=="Age") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_age <- ggplot(age, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.7 Probability of exclusion as a function of GAI

gai <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_gai <- ggplot(gai, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

p_legend = cowplot::get_legend(p_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

2b. Density plots of covariates used to fit GAMs (across included & excluded children)

Define theme for density plots of covariates across included and excluded children

2b.1 ADOS density

ddata <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data) %>% 
  filter(PrimaryDiagnosis=="ASD")

ddata$PrimaryDiagnosis <- droplevels(ddata$PrimaryDiagnosis)

d_ados=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .09), breaks=seq(0, .08, by=.02))+  
  labs(x='', y='Density')+
  scale_fill_manual(values = c("#FDE599"))+ 
  scale_color_manual(values = c("#E9D38D"))+ 
  den_theme

2b.2 SRS density

ddata <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_srs=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(0,max(srs$range)),breaks = seq(0, 100 , by = 50))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.3 Inattention density

ddata <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_in=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  labs(x='')+  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  theme(axis.title.y = element_blank())+
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.4 Hyperactivity/Impulsivity Density

ddata <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_hi=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.5 Motor Overflow Density

ddata <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_mo=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  theme(axis.title.y = element_blank())+
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.6 Age Density

ddata <- nested_gams %>% 
  filter(variable=="Age") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_age=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(8,13), breaks = seq(8, 13 , by = 1))+  
  scale_y_continuous(expand = c(0, 0), limits=c(0,.29), breaks=seq(0, .25, by = .05))+ 
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.7 GAI Density

ddata <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_gai=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .035), breaks=seq(0., .03, by=.01))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  labs(x='')+  
  theme(axis.title.y = element_blank())

hist_legend = cowplot::get_legend(d_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

combine gam plots with densities & print

top_row <- cowplot::plot_grid(p_ados, p_srs, p_in, p_hi, p_mo, p_age, p_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
bottom_row <- cowplot::plot_grid(d_ados, d_srs, d_in, d_hi, d_mo, d_age, d_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
## Warning: Removed 91 rows containing non-finite values (stat_density).
## Warning: Removed 19 rows containing non-finite values (stat_density).
png("Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))
dev.off()
## quartz_off_screen 
##                 2
#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))

#dev.off()

NOTE: 19 rows = # of participants with missing Motor Overflow scores (imputed for deconfounded group difference)

3. Impact of motion QC on phenotypic representation within diagnostic groups

3a. Split violin plots

My_Theme = theme_light()+theme(
  legend.title = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  axis.text.x = element_text(size = 5),
  axis.text.y = element_text(size = 8),
  strip.text.x = element_text(size = 12, face = "bold", color="black"),
  strip.text.y = element_text(size = 10, color="black"),
  strip.background = element_rect(fill="white"),
  plot.title = element_text(size = 9, hjust = 0.5))

3a.1 ADOS (ASD only) split violin

NOTE: Missing ADOS scores for one participant evaluated at CARD

ados <- aim1 %>% 
  filter(PrimaryDiagnosis=="ASD" & variable=="ADOS") %>% 
  dplyr::select(-c("PrimaryDiagnosis", "variable"))

ados_violin <- ggplot(ados, aes(Motion.Exclusion.Level, value, fill = Included, color = Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 1.5) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.text = element_text(size = 8))+
  theme(legend.position = "none")+
  ggtitle("ADOS")

3a.2 SRS split violin

srs <- filter(aim1G, variable=="SRS")
aim1_srs <- ggplot(srs, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~., scales = "fixed")+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  # geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("SRS")
  
v_legend = cowplot::get_legend(aim1_srs + guides(color = guide_legend(nrow = 2))+
                                 theme(legend.position = "left", 
                                       legend.text = element_text(size = 10), 
                                       legend.key.size=unit(.1, "in")))
## Warning: Removed 273 rows containing non-finite values (stat_ydensity).
## Warning: Removed 273 rows containing non-finite values (stat_summary).

## Warning: Removed 273 rows containing non-finite values (stat_summary).

NOTE: 273/3 = 91, # of participants missing SRS

with(dat3, table(PrimaryDiagnosis, is.na(SRS.Score)))
##                 
## PrimaryDiagnosis FALSE TRUE
##              TD    263   90
##              ASD   150    1

3a.3 Inattention Split Violin

inat <- filter(aim1G, variable=="Inattention")
paper_inat <- ggplot(inat, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Inattention")

3a.4 Hyperactivity/Impulsivity Spilt Violin

hyp <- filter(aim1G, variable=="Hyperactivity")
paper_hyp <- ggplot(hyp, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Hyperactivity")

3a.5 Motor Overflow Split Violin

overflow <- filter(aim1G, variable=="Motor Overflow")
aim1_of <- ggplot(overflow, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Motor Overflow")

3a.6 Age Split Violin

age <- filter(aim1G, variable=="Age")
aim1_age <- ggplot(age, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Age")

3a.7 GAI Split Violin

gai <- filter(aim1G, variable=="GAI")
aim1_gai <- ggplot(gai, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.position = "none")+
  ggtitle("GAI")

Combine split violins into one figure & save

3b. Mann-Whitney U tests to compare included vs excluded participants using lenient motion QC (13 tests)

#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Lenient"))

aim1MW <- aim1tib %>% 
  group_by(PrimaryDiagnosis, variable) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms. NOTE: ADOS only collected in ASD (9 tests)
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

#hypothesis: included children will be older and have greater GAI (4 tests)
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]
## # A tibble: 13 × 7
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable       statistic  p.value method  alternative  p.fdr
##    <fct>            <fct>              <dbl>    <dbl> <chr>   <chr>        <dbl>
##  1 ASD              ADOS               1780. 0.00921  Wilcox… less        0.0264
##  2 ASD              SRS                1720  0.00579  Wilcox… less        0.0264
##  3 ASD              Motor Overflow     1260  0.000948 Wilcox… less        0.0123
##  4 ASD              Inattention        2442. 0.642    Wilcox… less        0.642 
##  5 ASD              Hyperactivity      2419  0.606    Wilcox… less        0.642 
##  6 ASD              Age                2848  0.0216   Wilcox… greater     0.0352
##  7 ASD              GAI                2960  0.00655  Wilcox… greater     0.0264
##  8 TD               SRS                4562. 0.431    Wilcox… less        0.509 
##  9 TD               Motor Overflow     7086. 0.0569   Wilcox… less        0.0822
## 10 TD               Inattention        8004  0.268    Wilcox… less        0.349 
## 11 TD               Hyperactivity      7018. 0.0198   Wilcox… less        0.0352
## 12 TD               Age               10074. 0.0102   Wilcox… greater     0.0264
## 13 TD               GAI                9882  0.0202   Wilcox… greater     0.0352

3c. Mann-Whitney U tests to compare included vs excluded participants using strict motion QC (13 tests)

#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Strict"))

aim1MW <- aim1tib %>% 
  group_by(variable, PrimaryDiagnosis) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

#hypothesis: included children will be older and have greater GAI
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
  coefs = map(mwm, tidy)) %>% 
  unnest(coefs)

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]
## # A tibble: 13 × 7
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable       statistic p.value method   alternative  p.fdr
##    <fct>            <fct>              <dbl>   <dbl> <chr>    <chr>        <dbl>
##  1 ASD              ADOS               1350. 0.0236  Wilcoxo… less        0.0767
##  2 ASD              SRS                1312. 0.0176  Wilcoxo… less        0.0763
##  3 ASD              Motor Overflow     1132. 0.0437  Wilcoxo… less        0.0847
##  4 ASD              Inattention        1582  0.189   Wilcoxo… less        0.223 
##  5 ASD              Hyperactivity      1671  0.322   Wilcoxo… less        0.322 
##  6 ASD              Age                2221  0.0165  Wilcoxo… greater     0.0763
##  7 ASD              GAI                1893  0.280   Wilcoxo… greater     0.303 
##  8 TD               SRS                7274. 0.0456  Wilcoxo… less        0.0847
##  9 TD               Motor Overflow    13570. 0.149   Wilcoxo… less        0.194 
## 10 TD               Inattention       14006. 0.147   Wilcoxo… less        0.194 
## 11 TD               Hyperactivity     12199  0.00122 Wilcoxo… less        0.0159
## 12 TD               Age               16656. 0.0374  Wilcoxo… greater     0.0847
## 13 TD               GAI               16378. 0.0685  Wilcoxo… greater     0.111

4. Edgewise functional connectivity as a function of the covariates

Combine partial correlations with melted rs-fMRI usability and covariate info

startEdgeidx = which(names(dat3)=='r.ic1.ic2')
endEdgeidx = which(names(dat3)=='r.ic29.ic30')
names(dat3)[startEdgeidx:endEdgeidx]
##   [1] "r.ic1.ic2"   "r.ic1.ic4"   "r.ic1.ic8"   "r.ic1.ic13"  "r.ic1.ic14" 
##   [6] "r.ic1.ic15"  "r.ic1.ic17"  "r.ic1.ic19"  "r.ic1.ic21"  "r.ic1.ic22" 
##  [11] "r.ic1.ic24"  "r.ic1.ic25"  "r.ic1.ic26"  "r.ic1.ic27"  "r.ic1.ic28" 
##  [16] "r.ic1.ic29"  "r.ic1.ic30"  "r.ic2.ic4"   "r.ic2.ic8"   "r.ic2.ic13" 
##  [21] "r.ic2.ic14"  "r.ic2.ic15"  "r.ic2.ic17"  "r.ic2.ic19"  "r.ic2.ic21" 
##  [26] "r.ic2.ic22"  "r.ic2.ic24"  "r.ic2.ic25"  "r.ic2.ic26"  "r.ic2.ic27" 
##  [31] "r.ic2.ic28"  "r.ic2.ic29"  "r.ic2.ic30"  "r.ic4.ic8"   "r.ic4.ic13" 
##  [36] "r.ic4.ic14"  "r.ic4.ic15"  "r.ic4.ic17"  "r.ic4.ic19"  "r.ic4.ic21" 
##  [41] "r.ic4.ic22"  "r.ic4.ic24"  "r.ic4.ic25"  "r.ic4.ic26"  "r.ic4.ic27" 
##  [46] "r.ic4.ic28"  "r.ic4.ic29"  "r.ic4.ic30"  "r.ic8.ic13"  "r.ic8.ic14" 
##  [51] "r.ic8.ic15"  "r.ic8.ic17"  "r.ic8.ic19"  "r.ic8.ic21"  "r.ic8.ic22" 
##  [56] "r.ic8.ic24"  "r.ic8.ic25"  "r.ic8.ic26"  "r.ic8.ic27"  "r.ic8.ic28" 
##  [61] "r.ic8.ic29"  "r.ic8.ic30"  "r.ic13.ic14" "r.ic13.ic15" "r.ic13.ic17"
##  [66] "r.ic13.ic19" "r.ic13.ic21" "r.ic13.ic22" "r.ic13.ic24" "r.ic13.ic25"
##  [71] "r.ic13.ic26" "r.ic13.ic27" "r.ic13.ic28" "r.ic13.ic29" "r.ic13.ic30"
##  [76] "r.ic14.ic15" "r.ic14.ic17" "r.ic14.ic19" "r.ic14.ic21" "r.ic14.ic22"
##  [81] "r.ic14.ic24" "r.ic14.ic25" "r.ic14.ic26" "r.ic14.ic27" "r.ic14.ic28"
##  [86] "r.ic14.ic29" "r.ic14.ic30" "r.ic15.ic17" "r.ic15.ic19" "r.ic15.ic21"
##  [91] "r.ic15.ic22" "r.ic15.ic24" "r.ic15.ic25" "r.ic15.ic26" "r.ic15.ic27"
##  [96] "r.ic15.ic28" "r.ic15.ic29" "r.ic15.ic30" "r.ic17.ic19" "r.ic17.ic21"
## [101] "r.ic17.ic22" "r.ic17.ic24" "r.ic17.ic25" "r.ic17.ic26" "r.ic17.ic27"
## [106] "r.ic17.ic28" "r.ic17.ic29" "r.ic17.ic30" "r.ic19.ic21" "r.ic19.ic22"
## [111] "r.ic19.ic24" "r.ic19.ic25" "r.ic19.ic26" "r.ic19.ic27" "r.ic19.ic28"
## [116] "r.ic19.ic29" "r.ic19.ic30" "r.ic21.ic22" "r.ic21.ic24" "r.ic21.ic25"
## [121] "r.ic21.ic26" "r.ic21.ic27" "r.ic21.ic28" "r.ic21.ic29" "r.ic21.ic30"
## [126] "r.ic22.ic24" "r.ic22.ic25" "r.ic22.ic26" "r.ic22.ic27" "r.ic22.ic28"
## [131] "r.ic22.ic29" "r.ic22.ic30" "r.ic24.ic25" "r.ic24.ic26" "r.ic24.ic27"
## [136] "r.ic24.ic28" "r.ic24.ic29" "r.ic24.ic30" "r.ic25.ic26" "r.ic25.ic27"
## [141] "r.ic25.ic28" "r.ic25.ic29" "r.ic25.ic30" "r.ic26.ic27" "r.ic26.ic28"
## [146] "r.ic26.ic29" "r.ic26.ic30" "r.ic27.ic28" "r.ic27.ic29" "r.ic27.ic30"
## [151] "r.ic28.ic29" "r.ic28.ic30" "r.ic29.ic30"
signalFC <- dat3[, c(1,startEdgeidx:endEdgeidx)]

dat2 <- merge(aim1, signalFC, all=TRUE, by = "ID")
fcMelt <- reshape2::melt(dat2[,1:160],
                         id.vars=names(dat2)[1:7],
                         variable.name = "edge",
                         value.name = "fc")

4a. Run nested gams

fctib <- tibble(filter(fcMelt, Motion.Exclusion.Level!="None"))
fctib$Motion.Exclusion.Level <- droplevels(fctib$Motion.Exclusion.Level)

fcNest <- fctib %>% 
  filter(Included=="Included") %>% 
  group_by(variable, Motion.Exclusion.Level, edge) %>% 
  tidyr::nest()

#nested models
fc_gams <- fcNest %>% 
  mutate(model = map(data, ~mgcv::gam(fc~s(value), data = na.omit(.x), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE)) %>% 
  unnest(coefs)

4b. Plot histograms of edgewise p-values from GAMs

Histogram of edgewise p-values for partial correlations as a function of ADOS.

NOTE: TD scores = 0

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="ADOS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_ados_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  #scale_fill_manual(labels=c('TD','ASD'), values = c("#FDE599", "#9FB0CC"))+
  #scale_color_manual(labels=c('TD','ASD'), values = c("#E9D38D", "#8C9AB4"))+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  labs(x='p value', y='Count')+
  theme(axis.title.y = element_text(size = 9))+
  theme(axis.title.x = element_text(size = 7))

4b.1 Histogram of edgewise p-values for partial correlations as a function of SRS

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="SRS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_srs_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.2 Histogram of p-values for partial correlations as a function of inattention

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Inattention") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_in_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.3 Histogram of p-values for partial correlations as a function of hyperactivity

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Hyperactivity") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_hi_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Hyperactivity")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.4 Histogram of p-values for partial correlations as a function of motor overflow

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Motor Overflow") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_mo_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.5 Histogram of p-values for partial correlations as a function of age

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Age") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_age_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.6 Histogram of p-values for partial correlations as a function of GAI

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="GAI") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_gai_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

hist_p_legend <- cowplot::get_legend(hist_gai_p + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in"), legend.title = element_blank()))

4b.7 Combine histograms and print

fc_hist <- cowplot::plot_grid(hist_ados_p, hist_srs_p, hist_in_p, hist_hi_p, hist_mo_p, hist_age_p, hist_gai_p, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))

#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_hist_fc.png",width=8,height=3,units="in",res=200)
png("Figures/fig_hist_rfc_facet.png",width=8,height=3,units="in",res=200)

cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))
dev.off()
## quartz_off_screen 
##                 2
cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))